% SimGreedModelOnceFinal.m  - One iteration of the numerical Greed model to replicate Figure 1
% Guardado and Pennings (2023) "Shock Persistence and the Study of Armed Conflict: Empirical Biases and Some Remedies" International Interactions
% Contact: Steven Pennings (spennings@worldbank.org; steven.pennings@gmail.com)
%       Replicates Figure1A.pdf (Figure 1, Panel A) and Figure1BC.pdf (Figure 1, Panels B and C)
%       Calls solveGreedSS.m; GreedAR1.mod; GreedSV.mod
%       Uses Matlab R2018a (Optimization Toolbox; Statistics & Machine Learning toolbox); Dynare 5.3 (www.dynare.org)



%% Model set up %%
clear
clc
cd C:\Users\wb464833\Dropbox\SeasonalViolence\TheoryNote\Replication\GreedModel
addpath C:\dynare\5.3\matlab
 
% Parameters (we add the "1" as use this to keep variable later on")
global beta thetabar alpha phi psi  rho lambda 
beta=0.99;                  % Quarterly discount rate (in the manuscript we call this delta, because beta is the coefficient in a regression)
beta1=beta;
thetabar=1;                 % Steady State TFP level(not logged!)
thetabar1=thetabar;
alpha=0.5;                  % Share of labor in production.
alpha1=alpha;
phi=3;                      % Opportunity cost is -phi (phi>1 to ensure a unique equilibrium)
phi1=phi;
psi=0.0065;                 % Scaling factor on the probability of winning the war
psi1=psi;
lambda=1;                   % Efficiency cost of government run by the rebels (1 is default; shouldn't be too low>0.5 or 0.6)
lambda1=lambda;
rho=0.966                   % Quarterly persistence of coffee prices (rho^4 annualized) - used in simulations in Figure 1A
%rho=0.983                  % Quarterly persistence of oil prices     (rho^4 annualized)
rho1=rho;

[vbar fval] = fsolve(@(x) solveGreedSS(x),[0.07])
vbarstar=vbar;                                       % Steady state value of violence from Appendix 1.3.1
vbarstar1=vbar; 
p=psi*vbar^(1-1/phi) ;                               % Probability of winning 
pprime=(1-1/phi)*psi*vbar^(-1/phi);                  % Probability of winning derivative
wbar=alpha*thetabar*(1-vbar)^(alpha-1);             % 1. Definition of wages=MPL
vbarW=lambda*thetabar/(1-beta);                     % 3. Value Fn of winning  (lambda-=1)
vbarL=(wbar*(1-vbar)+beta*p*vbarW)/(1-beta*(1-p));  % 4. Value Fn of losing (out of power)  (lambda-=1) 
vbar1=vbar;
vbarW1=vbarW;
vbarL1=vbarL;

rng(0,'twister') % This is the default Dynare seed

%% Figure 1A: One simulation of the model with AR(1) shocks of the same persistence as coffee prices  %%
clearvars -except *1
dynare GreedAR1.mod noclearall
[B,BINT] = regress(v,[ones(1000,1) w]);
B_AR1reg1=B(2)    % Measured opportunity cost of -2.2 from regression on 1000 quarters of simulated data (quoted in text)
%[B,BINT] = regress(v_thetashock,[ones(40,1) w_thetashock]);
%B_AR1regIRF1=B(2) % Measured opportunity cost of -2.2 from regression on 40 quarters of IRF
trueoppcost1=-phi1                                              %Store AR(1) output to use later in figure
vAR1=v(1:40);                                                   %Store AR(1) output to use later in figure
wAR1=w(1:40);                                                   %Store AR(1) output to use later in figure
piprizeAR1=(vbarW*EvW(1:40)-vbarL*EvL(1:40))/(vbarW-vbarL);

F1=figure (1)
title('A. Persistent shocks (5yr half-life, based on coffee price shocks)');
hold on
plot2Av=plot(vAR1*100);
plot2Aw=plot(wAR1*100);
plot2Ac=plot(phi1*piprizeAR1*100);
hold off 
axis([1 20 -0.3*100 0.11*100])
set(plot2Aw,'Color','black','LineWidth',2)
set(plot2Av,'Color','blue','LineWidth',2, 'Marker', 'x')
set(plot2Ac,'Color','red','LineWidth',2, 'Marker', 'o')
ax = gca;
ax.XGrid = 'off';
ax.YGrid = 'on';
legend({'Violence','Wages','\phi\timesPrize from fighting  [goes in error term as unobserved]'},'Location','SouthWest','FontSize',13)
legend('boxoff')
xlabel('Quarters')
ylabel('Deviation from Steady State (%)')
F1.Position
F1.Position(3)=1.5*560
ytickformat('percentage')
print(F1,'Figure1A.pdf','-dpdf','-bestfit')

%% Figure 1B: One simulation of the model with transient AR(1) shocks (rho=0)  %%
rho=0;
rho1=rho;
clearvars -except *1
dynare GreedAR1.mod noclearall
vAR2=v(1:10);
wAR2=w(1:10);
piprizeAR1_2=(vbarW*EvW(1:40)-vbarL*EvL(1:40))/(vbarW-vbarL);
[B,BINT] = regress(v,[ones(1000,1) w]);
B_transitent_reg1=B(2)            %Measured opportunity cost of -3 from regression using transient shocks (true value)

F2=figure (2)
subplot(1,2,1)
title('B. Transient Shocks (\rho=0)');
hold on
plot2Bv=plot(vAR2*100);
plot2Bw=plot(wAR2*100);
plot2Bc=plot(phi1*piprizeAR1_2*100);
set(plot2Bw,'Color','black','LineWidth',2)
set(plot2Bv,'Color','blue','LineWidth',2, 'Marker', 'x')
set(plot2Bc,'Color','red','LineWidth',2,'Marker', 'o')
hold off 
axis([1 10 -0.1*100 0.11*100])
ax = gca;
ax.XGrid = 'off';
ax.YGrid = 'on';
xlabel('Quarters')
yticks([-0.1*100 -0.05*100 0*100 0.05*100 0.1*100])
ylabel('Deviation from Steady State (%)')
ytickformat('percentage')
% We will then save the figure after the seasonal shock simulations

%% Figure 1C: One simulation of the model with Seasonal shocks %%%%
clearvars -except *1 *2
dynare GreedSV.mod noclearall
[B,BINT] = regress(v_thetashock,[ones(40,1) w_thetashock]);
B_SVregIRF1=B(2)            %Measured opportunity cost of -3 from regression using seasonal shocks (true value)
v_thetashock1=v_thetashock;
w_thetashock1=w_thetashock;
picombSV_thetashock1=(vbarW*EvW_thetashock-vbarL*EvL_thetashock)/(vbarW-vbarL);

subplot(1,2,2)
title('C. Seasonal Shocks');
hold on
plot3Cv=plot(v_thetashock1*100);
plot3Cw=plot(w_thetashock1*100);
plot3Cc=plot(phi1*picombSV_thetashock1*100);
hold off 
axis([1 10 -0.1*100 0.11*100])
set(plot3Cw,'Color','black','LineWidth',2)
set(plot3Cv,'Color','blue','LineWidth',2, 'Marker', 'x')
set(plot3Cc,'Color','red','LineWidth',2, 'Marker', 'o')
xlabel('Quarters')
ax = gca;
ax.XGrid = 'off';
ax.YGrid = 'on';
yticks([-0.1*100 -0.05*100 0*100 0.05*100 0.1*100])
ylabel('Deviation from Steady State (%)')
ytickformat('percentage')
F2.Position
F2.Position(3)=1.5*560
print(F2,'Figure1BC.pdf','-dpdf','-bestfit')
